// Guardando funcoes para programa de tomografia Roche rtom_gsl_fdf*
// VERBOSE FLAGS DISPONIVEIS:
// - VERBOSE_BLURR: Mostra string de blurr gaussiano
// - VERBOSE_SLOW: Mostra atuacao de funcoes X2_slow e dX2_slow


//#include <omp.h>

bool global_my_gsl_err = false;

class GSL_Integration;

void my_handler(const char * reason,const char * file,int line,int gsl_errno)
{
	cerr << "# - [" << gsl_errno << "] my error handler called because " << reason;
	// faco variavel global apontar para erro
	global_my_gsl_err = true;
}

template <> void Roche::NewI <gsl_vector> (gsl_vector &g)
{
  if(g.size != star.theta.size()){
    ostringstream erro;
    erro << "# - Roche::NewI(const gsl_vector &)\n"
	  << "# -  Operacao nao pode ser realizada, imagens de tamanhos diferentes\n"
	  << "# - star.I[" << star.I.size() << "] | New.I[" << g.size << "]" << endl;
          aborta(erro.str());
  }
  double  qq = 1./roche_par.q; // lobo da secundaria
  cout << "aqui!\n";
  for(int i = 0; i < g.size; i++){
    double rlob = ache_r(star.theta[i], star.phi[i]);
    double x = star.px[i];// *rlob;
    double y = star.py[i];// *rlob;
    double z = star.pz[i];// *rlob;
    double nr = sqrt(x*x+y*y+z*z);
    double n_r = sqrt((1-x)*(1-x)+(y*y)+(z*z));
    nr = nr*nr*nr;
    n_r = n_r*n_r*n_r;
    double nx = -x/nr+(qq*(1-x)/n_r)+(qq+1)*(x)-qq;
    double ny = -y*(1./nr+qq/n_r-(qq+1));
    double nz = -z*(1./nr+qq/n_r);
    double np = sqrt(nx*nx+ny*ny+nz*nz);
    //star.I[w] = pow(np/roche_par.g_pole, roche_par.obs_grav);
    star.I[i] = gsl_vector_get(&g,i)*pow(np/roche_par.g_pole, roche_par.obs_grav);
    //star.I[i]=newi[i]*pow(np/roche_par.g_pole, roche_par.obs_grav);
  }
}

  template <> void Roche::NewI <const vector<double> >(const vector<double> &newi)
{
  if(newi.size() != star.theta.size()) aborta("Nova mapa de tamanho diferente do esperado!");
  cout << "aqui!\n";
  double  qq = 1./roche_par.q; // lobo da secundaria
  for(int i = 0; i < newi.size(); i++){
    double rlob = ache_r(star.theta[i], star.phi[i]);
    double x = star.px[i];// *rlob;
    double y = star.py[i];// *rlob;
    double z = star.pz[i];// *rlob;
    double nr = sqrt(x*x+y*y+z*z);
    double n_r = sqrt((1-x)*(1-x)+(y*y)+(z*z));
    nr = nr*nr*nr;
    n_r = n_r*n_r*n_r;
    double nx = -x/nr+(qq*(1-x)/n_r)+(qq+1)*(x)-qq;
    double ny = -y*(1./nr+qq/n_r-(qq+1));
    double nz = -z*(1./nr+qq/n_r);
    double np = sqrt(nx*nx+ny*ny+nz*nz);
    //star.I[w] = pow(np/roche_par.g_pole, roche_par.obs_grav);

    star.I[i]=newi[i]*pow(np/roche_par.g_pole, roche_par.obs_grav);
  }

  //for(int i = 0; i < newi.size(); i++)
  //  star.I[i]=newi[i];
};


int definicoes_obscurecimento(const int law);
vector<vector<double> > org_par(vector<vector<double> > &in_par); //organiza os
void test_args(int, char**);
void show_help(void);
double objective(const gsl_vector *, void *);
void objective_div(const gsl_vector * , void *, gsl_vector *);
void ob_ob_div(const gsl_vector * , void * PARAMS, double * , gsl_vector *);
//bool terminator(GAGeneticAlgorithm &);

struct Dados {
  vector<double> fase;
  vector<double> fluxo;
  vector<int> min_x; //guarda de onde comessa o espectro de verdade flux(x) > 1e-4
  vector<int> max_x;
  vector<Histograma> dados;
  vector<Histograma> err;
};

enum function_list {Fdefault, Fslow, Fslow_2, Fslow_3};

class DataKeeper {
  //Definicoes para leitura de arquivo...
  int __IFDNP__; // = 3;
  int *__IFDP__;// = {6, 6, 7};
  int __IFDNCI__;// = 12;
  int __IFDB__ ; // = 2;
  //  const int IFD_SL = 3;
  //const int IFD_SP = 5;
  
  Dados _input_;
  vector<vector<double> > _par_;
  vector<vector<double> > _f2;
  vector<vector<int> > _pixmap_;
  function_list _fchoice;
  double _pfac_;
  double _X2_;
  double _S_;
  double _X2c_;
  double _Sc_;
  double _mval_;
  double __ZERO_MOD_FLUX__;
  double *_mlflux; //Model Line Flux
  int _NGBLURR;
  int __nfuncs;
  double _mfac;
  int _Np_; // Number of profiles for orbital smearing
  int _NITER; // Number of iterations
  Roche _estrela_;
  
  gsl_vector * _Default_Map;
	
  public:
    DataKeeper(int argc, char** argv);
    gsl_vector * Default_Map(void){return _Default_Map;};
    //float objective(GAGenome &);
	void Recalc_MLFLUX(void);
	void Recalc_Default_Map(const gsl_vector &,const double);
	void Init_Default_Map(void);
	
	/************
	 Diferent chi-square functions:
	  _slow - temptative introduction of instrumental smearing (OUTDATED - OUTUSE)
	  _slow_2 - working implementation of instrumental smearing
	  _slow_3 - temptative introduction of orbital smearing
	 ************/
    double inline Xaim(void){return _par_[0][1];};	
	void inline X2_reset(void){_X2_=-1;};
    double inline X2c(void){return _X2c_;};
    double inline X2(void){return _X2_;};
    double X2(const double);
    double X2(const gsl_vector &);
    double X2_slow(const gsl_vector &);
    double X2_slow_2(const gsl_vector &);
	double X2_slow_3(const gsl_vector &g);

	/************
	 chi-square derivatives:
	 ************/
	
	void dX2(const gsl_vector &, gsl_vector &);
    void dX2_slow(const gsl_vector &, gsl_vector &);
	void dX2_slow_3(const gsl_vector &, gsl_vector &);
	
	/************
	 Diferent Entropy functions:
	 ************/
	double inline S(void){return _S_;};
    double inline Sc(void){return _Sc_;}
    double S(const gsl_vector &);
	double S_gauss_slow(const gsl_vector &,const double);
	double S_gauss_fast(const gsl_vector &);

	/************
	 Entropy derivatives:
	 ************/
	
    void dS(const gsl_vector &,gsl_vector &);
	void dS_gauss_fast(const gsl_vector &,gsl_vector &);
	
	
	double inline mfac(void){
		return _mfac;
	};
	
	int inline nIter(void){return _NITER;};
	void inline nIter(int n){_NITER = n;};
	
    void inline mval(const gsl_vector &image)
    {
      double mm=0.0;
      for(int i = 0; i < image.size;i++)
	mm+=gsl_vector_get(&image,i);
      _mval_=mm/double(image.size);
    };
    double inline mval(void){return _mval_;};
	
    double inline pfac(void){return _pfac_;};

	double inline FWHM(void){
		return _par_[0][2];
	};
    void inline pfac(double f){_pfac_ = f;};
    int inline img_sz(void){return _estrela_.npix();};
    int inline data_sz(void){return _input_.fase.size();};
    surface inline image(void){return _estrela_.Star();};
    surface inline image(gsl_vector &g)
    {
      _estrela_.NewI(g);
      return _estrela_.Star();
    };
	
    void show_par(ostream& out);
	
    void inline nblurr(int val){_NGBLURR = val;};
    int inline nblurr(void){return _NGBLURR;};
	
    function_list fchoice(void){return _fchoice;};
    void fchoice(function_list val){_fchoice=val;};

    void inline write_res(string& root)
    {
      //vector<double> vfseq = vector<double>(101);
      //for (int i = 0 ; i <= 100; i++) vfseq[i] = i/100.0;
      vector<double>& vfases = _input_.fase;
      //if (mode == 1) vfases = vfseq;
      double &sn=_par_[0][3];//Sinal-Ruido
      for(int i = 0; i < vfases.size(); i++){
        // _input_.dados.col1() -> devolve referencia para vetor da coluna 1
        vector<double> lin_form(3);
        lin_form[0]=-_par_[1][4];
        lin_form[1]=_par_[1][4];
        lin_form[2]=_par_[1][5];
		  
		Histograma* Vec_MLProf = new Histograma [_Np_];
		  
		for (int nhist = 0; nhist < _Np_; ++nhist) {
			double dphi = _Np_ > 1 ? _par_[0][5]/(_Np_-1) * nhist - _par_[0][5]/2.0 : 0.0;
			  
			Vec_MLProf[nhist] = _estrela_.rad_prof(vfases[i] + dphi , _input_.dados[i].col1(), _pixmap_[i+nhist*vfases.size()]);
			  
		}
		  
		vector<vector<double> > line_mod(2); //,MLProf.size());
		line_mod[0].resize(Vec_MLProf[0].size());
		line_mod[1].resize(Vec_MLProf[0].size());
		  
		double norm_fac = 0.0; //_input_.fluxo[i]/_mlflux[i]; //Fator de normalizacao de fluxo (Fdata/Fmod)
		  
		for(int ii = 0; ii < Vec_MLProf[0].size(); ii++){
			line_mod[0][ii] = Vec_MLProf[0].le1(ii);
			line_mod[1][ii] = 0.0; //MLProf.le2(ii);
			  
			for (int kk = 0; kk < _Np_; ++kk) {
				  
				line_mod[1][ii] += Vec_MLProf[kk].le2(ii);
			}
			line_mod[1][ii] /= double(_Np_);
			  
			norm_fac+=line_mod[1][ii];
			  
		}
/*********************************************************************************************************		  
		for (int ii = 0; ii < line_mod[0].size(); ii++){
			line_mod[1][ii] *= _input_.fluxo[i]/norm_fac;
			  
		}
*********************************************************************************************************/
		
		if(_par_[1][6] > _par_[1][5]){
          Convolve ff(line_mod,_f2,line_mod[0].size(),_f2[0].size());
          ff.xmin(line_mod[0][0]);
          ff.xmax(line_mod[0][line_mod[0].size()-1]);
          norm_fac = _input_.fluxo[i]/_mlflux[i]; //ff.integrate1(); //_mlflux[i];
          //cout << i << " - aqui!\n";
           gsl_function F;
           F.function = &cnv_fnc;
           F.params = &ff;
           GSL_Integration itr_conv(F);
            //itr_conv.xmin(ff.xmin());
            //itr_conv.xmax(ff.xmax());
            itr_conv.xmin(line_mod[0][4]); //(ff.xmin());
            itr_conv.xmax(line_mod[0][line_mod[0].size()-4]); //(ff.xmax());
			norm_fac = 0;
			double avg_itrerr = 0.;
			for (int ii = 0; ii < line_mod[0].size(); ii++){
				ff.px(line_mod[0][ii]);
				line_mod[1][ii] = itr_conv.integrate_cquad()/sqrt(M_PI)/_par_[1][6];//_par_[1][6]/sqrt(M_PI)*norm_fac;
				norm_fac += line_mod[1][ii];
				avg_itrerr += itr_conv.integrationErr();
				itr_conv.clear();
					//cout << line_mod[0][ii]
					//<< " " << line_mod[1][ii] << endl;
			}
			if (avg_itrerr/line_mod[0].size() < 1.0)
				cout << " " << avg_itrerr/line_mod[0].size() << " ";
			ff.clean();
		}

		  //norm_fac /= double(Vec_MLProf[0].size());
		  
		delete [] Vec_MLProf;
		  
		norm_fac /= double(line_mod[0].size());

		for (int ii = 0; ii < line_mod[0].size(); ii++)
			line_mod[1][ii] *= _input_.fluxo[i]/_mlflux[i];// norm_fac;
				  
        //excluo os 2 primeiros e os 2 ultimos pontos do espectro devido a aplicacao do blurring!

		ostringstream sname;
        sname << root << ".dat_" << vfases[i];
        ofstream out(sname.str().c_str());
        cout << "# - Writing resulting spectrum in " << sname.str() << "...\n";
        for(int j = _input_.min_x[i]; j <= _input_.max_x[i]; j++){
          double flux = _input_.dados[i].le2(j);
          double err2 = _input_.err[i].le2(j)*sn;
          out << line_mod[0][j] << " " << line_mod[1][j] << " " 
               << _input_.dados[i].le1(j) << " " << flux << " " << err2 <<endl;
        }
      }
    };
};

void test_args(int argc, char** argv)
{
  if(argc < 2){
    show_help();
    exit(0);
  }
  if(!strcasecmp(argv[1],"-h")) {
    show_help();
    exit(0);
  }
}

void show_help(void)
{    
  cout << "\nModo de utilizacao do programa:\n"
	<< "\n\trtom_gafit config [img_def, [NITER]]\n\n"
	<< "Programa para realizar tomografia Roche pelo alargamento Dopler\n"
	<< "das linhas de absorcao/emissao de estrelas peenchendo o lobo de\n"
	<< "Roche. O arquivo de entrada 'config' é a saida de dgrid_rt. Para mais\n" 
	<< "informacoes rode 'dgrid_rt -h'.\n\n"
	<< "Tiago Ribeiro - LinuxC++ - 2008\n\n"
	<< "Last revision 01/march/2012.";
}

double objective(const gsl_vector * g,void * PARAMS)
{
  DataKeeper * dkp = (DataKeeper *) PARAMS;

//  DataKeeper &dkp = *(DataKeeper *) gen.userData();
  double X2=0.0; // = dkp->X2(*g);
  switch (dkp->fchoice()){
    case Fdefault:
      X2=dkp->X2(*g);
      break;
    case Fslow:
      X2=dkp->X2_slow(*g);
      break;
	case Fslow_2:
		  X2=dkp->X2_slow_2(*g);
		  break;
	case Fslow_3:
		  X2=dkp->X2_slow_3(*g);
		  break;
		  
	default:
      aborta("Funcao X2 nao especificada!");
  }
  double Xaim = dkp->Xaim();
  //double S = dkp->S(*g);//dkp->S(*g);
  double S = dkp->S_gauss_fast(*g);//dkp->S(*g);
  
  //cout << "# - INTERNAR [S = " << S << "]" << endl;
	
  if (Xaim < 1) Xaim+=1.0;
  //Xaim+=1.0;
	double c1 = (X2-Xaim)/Xaim;
  //cout << S*dkp->pfac() << " " << X2 << " " << c1 << " " << c1-S*dkp->pfac() << endl;
	return c1-S*dkp->pfac();//dkp->pfac()*(c1*c1)-S;
  //return -S*dkp->pfac();//dkp->pfac()*(c1*c1)-S;
}

void objective_div(const gsl_vector * x, void * PARAMS, gsl_vector * div)
{
  DataKeeper * dkp = (DataKeeper *) PARAMS;
  gsl_vector* divX2 = gsl_vector_alloc((*x).size);
  
  switch (dkp->fchoice()){
    case Fdefault:
      dkp->dX2(*x,*divX2);
      break;
    case Fslow:
      dkp->dX2_slow(*x,*divX2);
      break;
	  case Fslow_2:
		  dkp->dX2_slow(*x,*divX2);
		  break;
	  case Fslow_3:
		  dkp->dX2_slow_3(*x,*divX2);
		  break;
	  default:
      aborta("Funcao dX2 nao especificada!");
  }
  
  gsl_vector* divS = gsl_vector_alloc((*x).size);
  dkp->dS_gauss_fast(*x,*divS);
  //dkp->dS(*x,*divS);
  
  for(int i = 0; i < (*x).size; i++){
	  gsl_vector_set(div,i,-gsl_vector_get(divX2,i)+gsl_vector_get(divS,i)*dkp->pfac());
	  //cout << -gsl_vector_get(divX2,i)+gsl_vector_get(divS,i)*dkp->pfac() << " " 
	  //<< -gsl_vector_get(divX2,i) << " " << gsl_vector_get(divS,i)*dkp->pfac() << endl;
  }
	  //gsl_vector_set(div,i,gsl_vector_get(divS,i));
      
	//gsl_vector_set(div,i,gsl_vector_get(divS,i)*dkp->pfac());
    //gsl_vector_set(div,i,0.0);
  gsl_vector_free(divX2);
  gsl_vector_free(divS);
}

void ob_ob_div(const gsl_vector * x, void * PARAMS, double *ob , gsl_vector *div)
{
//  DataKeeper & dkp = (DataKeeper &) PARAMS;
  *ob=objective(x,PARAMS);
  objective_div(x,PARAMS,div);
}


DataKeeper::DataKeeper(int argc,char** argv)
{
  if(argc < 2) aborta("DataKeeper - Numero de parametros errado...");
  
  __IFDNP__ = 3;
  __IFDP__ = new int [__IFDNP__];
  __IFDP__[0] = 6;
  __IFDP__[1] = 7;
  __IFDP__[2] = 7;
  __IFDNCI__ = 12;
  __IFDB__ = 2;
  _NGBLURR = 200;
  __nfuncs = 2; //numero de funcoes
  _fchoice = Fdefault;
  __ZERO_MOD_FLUX__ = -1e-10;
  _X2_=-1.0;
  _S_=1.0;
  _mfac = ( argc > 2 ) ? atof(argv[2]) : 1.0;
  _NITER = ( argc > 3 ) ? atof(argv[3]) : 99;
    
  DataFileReader read(argv[1]);
  
  _par_.resize(__IFDNP__+1);
  
  for(int i = 0; i < __IFDNP__; i++){
    ++read;
    if(read.numberOfColumns() == __IFDP__[i]){
      double data_in;
      for(int j = 0; j < __IFDP__[i]; j++){
	 data_in = read.readColumn<double>(j+1);
	if(i == 1 && j == 0) data_in *= M_PI/180;
	_par_[i].push_back(data_in);
      }
    }
    else { 
      ostringstream erro;
      erro << "Numero de parametros errado\n\tArquivo: " 
	  << argv[1] << "\n\tLinha  : " << read.lineNumber();
      aborta(erro.str());
    }
/*
    if(read.eof()) {
      ostringstream erro;
      erro << "Numero de parametros errado\n\tArquivo: " 
	  << argv[1] << "\n\tLinha  : " << read.lineNumber();
      aborta(erro.str());
    }
    */
  }
  
    //Lendo coeficientes do obscurecimento de borda
  ++read;
  
  if(read.numberOfColumns()==definicoes_obscurecimento(int(_par_[2][__IFDB__]))){
    double data_in;
    for(int i = 0; i < definicoes_obscurecimento(int(_par_[2][__IFDB__])); i++){
       data_in = read.readColumn<double>(i+1);
      _par_[3].push_back(data_in);
    }
  }
  else {
    ostringstream erro;
    erro << "Numero de parametros errado\n\tArquivo: " 
	<< argv[1] << "\n\tLinha  : " << read.lineNumber();
    aborta(erro.str());
  }

  ++read;
  
  int LineFlag = 0;
  if(read.lineNumber() > 0)
    LineFlag = read.lineNumber();

  for(int i = 0; i < read.numberOfColumns(); i++){ // le lista ate o fim
    string filename=read.readColumn<string>(i+1);
    //fin >> filename;
    DataFileReader lpf(filename.c_str()); // abre arquivo da lista!

    //dar um geito de pegar fase do nome do arquivo e...

    string::size_type las = filename.find_last_of("_");
    string Sfase (filename, las+1);
    double fase = atof(Sfase.c_str());

    //cout << "# " << las << " - " << Sfase << " = " << fase << endl;

    _input_.fase.push_back(fase);

    vector<vector<double> > rdata(2);
    vector<vector<double> > err_rdata(2);
    while(!++lpf){ // le arquivo aberto ate o fim
      if(lpf.numberOfColumns()!= 3){ 
	ostringstream erro; 
	erro << "Numero de parametros errado\n\tArquivo: " 
	    << filename << "\n\tLinha  : " << lpf.lineNumber();
	aborta(erro.str());
      }
//      int VL = read.lineNumber()-LineFlag;
      rdata[0].push_back(lpf.readColumn<double>(1)*1e5);
      rdata[1].push_back(lpf.readColumn<double>(2) );
      err_rdata[0].push_back(lpf.readColumn<double>(1)*1e5);
      err_rdata[1].push_back(lpf.readColumn<double>(3));
    }
    Histograma rhist(rdata);
    _input_.dados.push_back(rhist);
    Histograma err_rhist(err_rdata);
    _input_.err.push_back(err_rdata);
    lpf.close();
    //++fin;
  }
  ++read;
  //Lendo fluxo dos espectros
  _input_.fluxo.resize(read.numberOfColumns());

  for(int i = 0; i < read.numberOfColumns(); i++)
    _input_.fluxo[i] = atof(read.readColumn<string>(i+1).c_str());
	
  //++read;
  _estrela_ = Roche(int(_par_[0][0]),org_par(_par_));
  //taking out orbital velocity from the line
  _estrela_.vorb(true);
  _estrela_.refaz_super();
	_estrela_.CalcPixSize();
  cout<<"# - Npix = "<<_estrela_.npix() << endl 
	<< "# - Apix = " << _estrela_.pix_size() <<endl;
  LineFlag=0;


  //Calcula numero de perfis de linha nescessarios
	
  _Np_ = ceil(_par_[0][5] / ( _par_[1][5]/1e5/(2.*M_PI*_estrela_.K2()) ) ); //Numero de perfis para combinar
	
  if (_Np_ < 2 ) { //This is as lower as it can go!
	  //Changing fchoice
	  //fchoice(Fslow_2);
	  _Np_ = 1;

  }
	
  cout << _Np_ << endl;

  /************
   Reservando espaco para mapeamento da posicao dos pixels no histograma. O mapa possui o mesmo numero de elementos
   das fases (cada pixel em fase) vezes o numero de perfis a serem somados. Cada mapa eh subsequente ao proximo de modo que
   eu leio todo o mapa 1 primeiro depois leio todo o mapa 2 e assim sucessivamente.
   ************/
  _pixmap_.resize(_input_.fase.size() * _Np_); 
	
  while(!++read){
    //nao vou ler a primeira coluna por que ela contem apenas o indice do pixel...
 
    if(read.numberOfColumns()-1 !=  _input_.fase.size() * _Np_){
		ostringstream erro; 
		erro << "\n\tNumero de pixels em mapa(" << read.numberOfColumns()-1 << ") diferente do esperado ("
		<<_input_.fase.size() * _Np_<<")...\n\t"
		<< "Verifique se texp/Porb em dgrid_rt eh igual ao utilizado para rodar rtom!\n";
		aborta(erro.str());
		
	}

    for(int i = 1; i < read.numberOfColumns(); i++)
      _pixmap_.at(i-1).push_back(read.readColumn<int>(i+1));
    ++LineFlag;
  }

  read.close();
  _f2.resize(2);
  //nblurr(200);
  //cria gaussiana com nblurr() pontos (default = 200)
  _f2[0].resize(int(floor(_par_[1][4]*nblurr()/_par_[1][5])));
  _f2[1].resize(int(floor(_par_[1][4]*nblurr()/_par_[1][5])));
  //double dxx = (in_par[0][5]*2.0/in_par[0][6]);
  //cout << dxx << endl;
  int itr=0;
  for ( double x0 = -_par_[1][4]; x0 < _par_[1][4]; x0+=_par_[1][5]/nblurr()*2.0){
    _f2[0][itr] = x0;
    _f2[1][itr] = exp(-(x0/_par_[1][6])*(x0/_par_[1][6]));
#ifdef VERBOSE_BLURR
    cout << _f2[0][itr] << " " << _f2[1][itr] << endl;
#endif
    itr++;
    if(itr >= _f2[0].size()) break;
  }

  //exit(0);
  //Pegando fluxo do modelo para utilizar como template
  _mlflux = new double [_input_.fluxo.size()];
  for( int i = 0; i < _input_.fluxo.size(); i++){
    Histograma MLProf = _estrela_.rad_prof(_input_.fase[i], _input_.dados[i].col1(), _pixmap_[i]);
    vector<vector<double> > mlp_vec(2);
    mlp_vec[0] = MLProf.col1();
    mlp_vec[1] = MLProf.col2();

    //GSL_Interp<vector<vector<double> > > itr_data(mlp_vec,MLProf.size());
    //_mlflux[i] = itr_data.integrate();
    double avg=0.;
    for (int jj = 0 ; jj < mlp_vec[0].size(); jj++)
		avg+=mlp_vec[1][jj];
    avg/=(double)mlp_vec[0].size();
	_mlflux[i] = avg;
//    int min = 0;
//    int max = 0;
	_input_.min_x.push_back(1);
	_input_.max_x.push_back(mlp_vec[0].size()-2);

	  /*
    for (int jj = 0 ; jj < mlp_vec[0].size(); jj++){
		
      if(mlp_vec[1][jj] > avg*1e-6 && min == 0) {
        _input_.min_x.push_back(jj-2);
        min = 1;

      }
      if(mlp_vec[1][mlp_vec[0].size()-1-jj] > avg*1e-6 && max == 0) {
        _input_.max_x.push_back(mlp_vec[0].size()-1-jj+2);
        max = 1;

      }
      if(min == 1 && max ==1)break;
		
    }
	 */
  }
/*
  for( int i = 0; i < _input_.fluxo.size(); i++){
    //string input = lis_in2.readColumn<string>(1);
    //saida << input << " ";
    //Calcular fluxo do espectro
    //DataFileReader data(input);
    vector<vector<double> > dados;
    dados.push_back(_input_.dados[i].col1());
    dados.push_back(_input_.dados[i].col2());
    //while(!++data){
    //  dados[0].push_back(data.readColumn<double>(1));
    //  dados[1].push_back(data.readColumn<double>(2));
    //}

	  GSL_Interp<vector<vector<double> > > itr_data(dados,dados[0].size());
    itr_data.xmin(_input_.dados[i].le1(_input_.min_x[i])); //(line_mod[0][4]); //(ff.xmin());
    itr_data.xmax(_input_.dados[i].le1(_input_.max_x[i])); //(line_mod[0][line_mod[0].size()-4]); //(ff.xmax());

    _input_.fluxo[i] = itr_data.integrate();
    //data.close();
    //lis_in2.seteol();
    //++lis_in2;
  }
 */
	Init_Default_Map();

}

void DataKeeper::Recalc_Default_Map(const gsl_vector &g,const double FWHM)
{
	surface Srf = image();
#pragma omp parallel for
	for(int i = 0; i < g.size; i++){
		double* pixlist;
		//cout << "# aqui 1\n";
		int npix = _estrela_.close_pix(Srf.px[i],Srf.py[i],Srf.pz[i],FWHM, pixlist);
		//cout << "# aqui 2\n";
		double Mj = 0.0;
		for( int j = 0; j < npix; j++)
			Mj += gsl_vector_get(&g,pixlist[j]);
		delete pixlist;
		//cout << "# aqui 3 - " << i << " / " << g.size << " / " << npix << "\n";
		Mj /= double(npix);
		gsl_vector_set(_Default_Map, i, Mj);
	}
}

void DataKeeper::Init_Default_Map(void)
{
	_Default_Map =	gsl_vector_alloc(img_sz());
}


void DataKeeper::Recalc_MLFLUX(void)
{

	for( int i = 0; i < _input_.fluxo.size(); i++){
		Histograma MLProf = _estrela_.rad_prof(_input_.fase[i], _input_.dados[i].col1(), _pixmap_[i]);
		vector<vector<double> > mlp_vec(2);
		mlp_vec[0] = MLProf.col1();
		mlp_vec[1] = MLProf.col2();
		//GSL_Interp<vector<vector<double> > > itr_data(mlp_vec,MLProf.size());
		//_mlflux[i] = itr_data.integrate();
		double avg=0;
		for (int jj = 0 ; jj < mlp_vec[0].size(); jj++)avg+=mlp_vec[1][jj];
		avg/=mlp_vec[0].size();
		_mlflux[i] = avg;
		int min = 0;
		int max = 0;
		/***********************************************
		for (int jj = 0 ; jj < mlp_vec[0].size(); jj++){
			if(mlp_vec[1][jj] > avg*1e-6 && min == 0) {
				_input_.min_x.push_back(jj-2);
				min = 1;
				
			}
			if(mlp_vec[1][mlp_vec[0].size()-1-jj] > avg*1e-6 && max == 0) {
				_input_.max_x.push_back(mlp_vec[0].size()-1-jj+2);
				max = 1;
				
			}
			if(min == 1 && max ==1)break;
		}
		***********************************************/
	}
	
}


double DataKeeper::X2(const gsl_vector &g)
{
	//Copia nova distribuicao de brilho para estrela...
	_estrela_.NewI(g);
#ifdef VERBOSE_FAST
	cout << "# - X2_fast started...\n";
#endif
	double X2 = 0;
	double NP = 0;
	double &sn=_par_[0][3];
	
	for(int i = 0; i < _input_.fase.size(); i++){
		double X2_phase = 0.0;
		// _input_.dados.col1() -> devolve referencia para vetor da coluna 1
		Histograma MLProf = _estrela_.rad_prof(_input_.fase[i], _input_.dados[i].col1(), _pixmap_[i]);
		
		//Preparando linha do modelo para calcular X2
		vector<vector<double> > line_mod(2);
		line_mod[0].resize(MLProf.size());
		line_mod[1].resize(MLProf.size());
		
		for(int ii = 0; ii < MLProf.size(); ii++){
			line_mod[0][ii] = MLProf.le1(ii);
			line_mod[1][ii] = MLProf.le2(ii);
		}
		
		double norm_fac = _input_.fluxo[i]/_mlflux[i];; //Fator de normalizacao de fluxo (Fdata/Fmod)
		
		for (int ii = 0; ii < line_mod[0].size(); ii++){
			line_mod[1][ii] *= norm_fac;
		}
		//Linha modelo preparada... Calculando X2
		
		int Nprof = 0; //conta o numero de pontos do perfil que entram na conta
		//for(int j = _input_.min_x[i]; j <= _input_.max_x[i]; j++){
		for(int j = 0; j <= _input_.dados[i].size(); j++){	
			double Vlow = MLProf.le1(j)-MLProf.step_down(j);
			double Vhi = MLProf.le1(j)+MLProf.step_up(j);
			
			double vel = _input_.dados[i].le1(j);
			
			if(Vlow <= vel && vel <= Vhi){
				double flux = _input_.dados[i].le2(j);
				double err2 = _input_.err[i].le2(j)*sn;
				//double err2 = _input_.dados[i].le2(j)*sn;
				err2*=err2;
				if(err2 == 0.0) err2=1.0;
				Nprof++;
				X2_phase +=(line_mod[1][j]-flux+_par_[2][6])*(line_mod[1][j]-flux+_par_[2][6])/err2;
				//X2+= (line_mod[1][j]-flux+_par_[2][6])*(line_mod[1][j]-flux+_par_[2][6])/err2;
#ifdef VERBOSE_FAST_ACTION
				cout << line_mod[0][j] << " " << line_mod[1][j] << " " 
				<< _input_.dados[i].le1(j) << " " << flux << " " << sqrt(err2) <<endl;
#endif
			}
/*			else if (MLProf.le2(j) > 0){
				cout << Vlow << "  <  " << vel << " < " << Vhi <<endl;
				cout << line_mod[0][j] << " " << line_mod[1][j] << " " 
				<< MLProf.le1(j) << " " << MLProf.le2(j) << endl;
				aborta("Problema no mapemento de pixels");
			}
*/			
		}
		
#ifdef VERBOSE_FAST_ACTION
		cout<<endl;
		exit(0);
#endif
		X2 += X2_phase/Nprof;
		NP+=Nprof; //double(_input_.dados[i].size());
		//cout << X2_phase/Nprof << endl;
		//exit(0);
	}
	X2/=_input_.fase.size();
	//exit(0);
	if((_X2_ < 0) ||(_X2_ > X2)) _X2_=X2; //guarga o menor valor de X2 em _X2_
	_X2c_=X2;
#ifdef VERBOSE_FAST
	cout << "# - X2_fast ended...\n";
#endif
	
	return X2;
	
	/* Primeira versao! obsoleta!
	 
  //Copia nova distribuicao de brilho para estrela...
  _estrela_.NewI(g);
  
  double X2 = 0;
  double NP = 0;
  double &sn=_par_[0][3];

  for(int i = 0; i < _input_.fase.size(); i++){
    // _input_.dados.col1() -> devolve referencia para vetor da coluna 1
    Histograma MLProf = _estrela_.rad_prof(_input_.fase[i], _input_.dados[i].col1(), _pixmap_[i]);
    
    for(int j = 0; j < MLProf.size(); j++){
      
      double Vlow = MLProf.le1(j)-MLProf.step_down(j);
      double Vhi = MLProf.le1(j)+MLProf.step_up(j);

      double vel = _input_.dados[i].le1(j);
      
      if(((Vlow <= vel && vel <= Vhi) && (MLProf.le2(j) > 0)&& (_input_.dados[i].le2(j) > 1e-10))){
	double flux = _input_.dados[i].le2(j);
	double err2 = _input_.err[i].le2(j); //flux/sn;
	err2*=err2;
	if(err2 == 0.0) err2=1.0;
	X2+= (MLProf.le2(j)-flux+_par_[2][6])*(MLProf.le2(j)-flux+_par_[2][6])/err2;
	//cout << _input_.dados[i].le1(j) << " "<< MLProf.le2(j) << " "
	//    << MLProf.le1(j) << " "<< flux << " " 
	//    << (MLProf.le2(j)-flux+_par_[2][6])*(MLProf.le2(j)-flux+_par_[2][6])/err2 <<endl;		
      }
      //else if(MLProf.le2(j) > 0) aborta("Problema no mapemento de pixels");
      
    }
//    cout << endl;

    NP+=double(_input_.dados[i].size());
  }
  X2/=NP;
  
  if((_X2_ < 0) ||(_X2_ > X2)) _X2_=X2; //menor
  _X2c_=X2; //current
  //cout<<X2<<endl;
  //exit(0);
  return X2;*/
}

double DataKeeper::X2(const double val)
{
  //Copia nova distribuicao de brilho para estrela...
//  GARealGenome & gen = (GARealGenome &)g;
  const vector<double> ni(img_sz(),val);

  _estrela_.NewI(ni);
  
  double X2 = 0;
  double NP = 0;
  double &sn=_par_[0][3];

  for(int i = 0; i < _input_.fase.size(); i++){
    // _input_.dados.col1() -> devolve referencia para vetor da coluna 1
    Histograma MLProf = _estrela_.rad_prof(_input_.fase[i], _input_.dados[i].col1(), _pixmap_[i]);
    
    for(int j = 0; j < MLProf.size(); j++){
      
      double Vlow = MLProf.le1(j)-MLProf.step_down(j);
      double Vhi = MLProf.le1(j)+MLProf.step_up(j);

      double vel = _input_.dados[i].le1(j);
      
      if(Vlow <= vel && vel <= Vhi){
	double flux = _input_.dados[i].le2(j);
	double err2 = flux/sn;
	err2*=err2;
	if(err2 == 0.0) err2=1.0;
	X2+= (MLProf.le2(j)-flux+_par_[2][6])*(MLProf.le2(j)-flux+_par_[2][6])/err2;

		
      }
      else if(MLProf.le2(j) > 0) aborta("Problema no mapemento de pixels");
      
    }

    NP+=double(_input_.dados[i].size());
  }
  X2/=NP;
  
  if((_X2_ < 0) ||(_X2_ > X2)) _X2_=X2;
  _X2c_=X2;
  return X2;
	
}

void DataKeeper::dX2(const gsl_vector &x, gsl_vector &rdX2)
{
	
	_estrela_.NewI(x);
#ifdef VERBOSE_FAST
	cout << "# - dX2_fast started...\n";
#endif
	
	//gsl_vector* rdX2 = gsl_vector_alloc(x.size); //vetor de retorno
	
	double &sn=_par_[0][3];//Sinal-Ruido
	double ABS_FLAG = _par_[0][4]; //Flag de absorcao
	//float ABS_FLAG = 1.0; // -1 p/ emissao e 1 para absorcao
	if (abs(ABS_FLAG) != 1.0) ABS_FLAG = ABS_FLAG/abs(ABS_FLAG);
	
	//Modeled Line Profile:
	vector<vector<double> > line_mod(_input_.fase.size());
	
	for(int k = 0; k < _input_.fase.size(); k++){
		Histograma MLProf = _estrela_.rad_prof(_input_.fase[k], _input_.dados[k].col1(), _pixmap_[k]);
		// aplicando blurring instrumental
		vector<double> velvec;
		velvec.resize(MLProf.size());
		line_mod[k].resize(MLProf.size());
		//cout << "#- Aqui 0\n";
		vector<vector<double> > line_mod2(2);
		line_mod2[0].resize(MLProf.size());
		line_mod2[1].resize(MLProf.size());
		
		for(int ii = 0; ii < MLProf.size(); ii++){
			velvec[ii] = MLProf.le1(ii);
			line_mod[k][ii] = MLProf.le2(ii);
			line_mod2[0][ii]=velvec[ii];
			line_mod2[1][ii]=line_mod[k][ii];
		}

		double norm_fac;

		norm_fac = _input_.fluxo[k]/_mlflux[k];
		
		for (int ii = 0; ii < velvec.size(); ii++){
			line_mod[k][ii] *= norm_fac;
		}
			
		
	}
	
	for(int i = 0; i < x.size; i++){
		
		double X2 = 0;
		for(int k = 0; k < _input_.fase.size(); k++){
			double vis = _estrela_.VFunc(i, _input_.fase[k]);
			if(vis > 0){ 
				// Posicao do pixel no histograma de velocidade 
				int HistPixPos = _pixmap_[k][i];
				//excluo os 2 primeiros e os 2 ultimos pontos do espectro devido a aplicacao do blurring!
				//if((2 <= HistPixPos) && (HistPixPos < _input_.dados.size()-2)){
				double flux = _input_.dados[k].le2(HistPixPos);
				double err2 = _input_.err[k].le2(HistPixPos)*sn; //flux/sn;
				//double err2 = _input_.dados[k].le2(HistPixPos)*sn; //flux/sn;
				//err2*=err2;
				if(err2 == 0.0) err2=1.0;
				X2+= 2.0*vis*(line_mod[k][HistPixPos]-flux-_par_[2][6]); //err2;
				//cout << line_mod[k][HistPixPos] << "  " <<flux << endl;
				//}	
			}
		}
		//cout << X2 << endl;
		//Para abcorcao tem que t
		
		gsl_vector_set(&rdX2, i, (ABS_FLAG*X2/double(_input_.fase.size())));
		
	}
#ifdef VERBOSE_FAST
	cout << "# - dX2_fast ended...\n";
#endif
	
	//return rdX2;
	
	
	/* Antigo e obseleto! 
  _estrela_.NewI(x);
  
  //gsl_vector* rdX2 = gsl_vector_alloc(x.size); //vetor de retorno
  
  double &sn=_par_[0][3];//Sinal-Ruido
  
  for(int i = 0; i < x.size; i++){

    double X2 = 0;
    for(int k = 0; k < _input_.fase.size(); k++){
      double vis = _estrela_.VFunc(i, _input_.fase[k]); 
      if(vis > 0){ 
	Histograma ModLineProf = _estrela_.rad_prof(_input_.fase[k], _input_.dados[k].col1(), _pixmap_[k]);
	// Posicao do pixel no histograma de velocidade 
	int HistPixPos = _pixmap_[k][i];
	
	if((0 <= HistPixPos) && (HistPixPos < _input_.dados.size())){
	  double flux = _input_.dados[k].le2(HistPixPos);
	  double err2 = _input_.err[k].le2(HistPixPos); //flux/sn;
	  err2*=err2;
	  if(err2 == 0.0) err2=1.0;
	  X2+= 2*vis*(ModLineProf.le2(HistPixPos)-flux+_par_[2][6])/err2;
	}	
      }
    }
    gsl_vector_set(&rdX2, i, (-X2/double(_input_.fase.size())));
  }
  //return rdX2;
	*/
}

//Slow functions apply instrumental blurring! Use it if instrumental bluring is considerable!
double DataKeeper::X2_slow(const gsl_vector &g)
{
  //Copia nova distribuicao de brilho para estrela...
  _estrela_.NewI(g);
#ifdef VERBOSE_SLOW
  cout << "# - X2_slow started...\n";
#endif
  double X2 = 0;
  double NP = 0;
  double &sn=_par_[0][3];

  for(int i = 0; i < _input_.fase.size(); i++){
	  double X2_phase = 0.0;
    // _input_.dados.col1() -> devolve referencia para vetor da coluna 1
    vector<double> lin_form(3);
    lin_form[0]=-_par_[1][4];
    lin_form[1]=_par_[1][4];
    lin_form[2]=_par_[1][5];
    Histograma MLProf = _estrela_.rad_prof(_input_.fase[i], _input_.dados[i].col1(), _pixmap_[i]);

    // aplicando blurring instrumental
    vector<vector<double> > line_mod(2);
    line_mod[0].resize(MLProf.size());
    line_mod[1].resize(MLProf.size());
    //cout << "#- Aqui 0\n";
    for(int ii = 0; ii < MLProf.size(); ii++){
      line_mod[0][ii] = MLProf.le1(ii);
      line_mod[1][ii] = MLProf.le2(ii);
      //cout << line_mod[0][ii] << "  " << line_mod[1][ii] << " " << _input_.fase[i] <<endl;
    }

    //cout << endl;
    //cout << "#- Aqui 1\n";
    //tmp_vec=vector<double> (lp.col2());
    //line_mod[1].push_back(tmp_vec);
    double norm_fac = 0.0; //Fator de normalizacao de fluxo (Fdata/Fmod)
    //se tiver blurring aplica...
    if(_par_[1][5] > 0){
#ifdef VERBOSE_SLOW
		cout << "# - X2_slow - Convolution started...\n";
#endif
		
      Convolve ff(line_mod,_f2,line_mod[0].size(),_f2[0].size());
      ff.xmin(_input_.dados[i].le1(_input_.min_x[i]));//line_mod[0][0]);
      ff.xmax(_input_.dados[i].le1(_input_.max_x[i]));//line_mod[0][line_mod[0].size()-1]);
      norm_fac = _input_.fluxo[i]/ff.integrate1(); //_mlflux[i]; /
      //cout << i << " - aqui!\n";
      for (int ii = 0; ii < line_mod[0].size(); ii++){
        gsl_function F;
        F.function = &cnv_fnc;
        F.params = &ff;
        GSL_Integration itr_conv(F);
        //itr_conv.xmin(ff.xmin());
        //itr_conv.xmax(ff.xmax());
        itr_conv.xmin(_input_.dados[i].le1(_input_.min_x[i])); //(line_mod[0][4]); //(ff.xmin());
        itr_conv.xmax(_input_.dados[i].le1(_input_.max_x[i])); //(line_mod[0][line_mod[0].size()-4]); //(ff.xmax());
        ff.px(line_mod[0][ii]);
        line_mod[1][ii] = itr_conv.integrate_qag()/_par_[1][5]/sqrt(M_PI)*norm_fac;
        itr_conv.clear();
        //cout << line_mod[0][ii]
        //<< " " << line_mod[1][ii] << endl;
#ifdef VERBOSE_SLOW
		  cout << "# - X2_slow - " << ii <<" x " <<line_mod[0].size()<< endl;
#endif
		  
		  
      }
      //norm_fac/=(double(line_mod[0].size())*ff.integrate1());
      //for (int ii = 0; ii < line_mod[0].size(); ii++) line_mod[1][ii]*norm_fac;
      ff.clean();
    }
	else{
		Convolve ff(line_mod,_f2,line_mod[0].size(),_f2[0].size());
		ff.xmin(_input_.dados[i].le1(_input_.min_x[i]));
		ff.xmax(_input_.dados[i].le1(_input_.max_x[i]));
		norm_fac = _input_.fluxo[i]/ff.integrate1(); //_mlflux[i];
		
		for (int ii = 0; ii < line_mod[0].size(); ii++){
			line_mod[1][ii] *= _input_.fluxo[i]/ff.integrate1();
		}
		
	}
	  
	  
#ifdef VERBOSE_SLOW
	  cout << "# - X2_slow - Convolution ended...\n";
#endif
	  
    //excluo os 2 primeiros e os 2 ultimos pontos do espectro devido a aplicacao do blurring!

	  int Nprof = 0; //conta o numero de pontos do perfil que entram na conta
    for(int j = _input_.min_x[i]; j <= _input_.max_x[i]; j++){
      
      double Vlow = MLProf.le1(j)-MLProf.step_down(j);
      double Vhi = MLProf.le1(j)+MLProf.step_up(j);

      double vel = _input_.dados[i].le1(j);
      
      if(Vlow <= vel && vel <= Vhi){
	double flux = _input_.dados[i].le2(j);
	double err2 = _input_.err[i].le2(j)*sn;
	//double err2 = _input_.dados[i].le2(j)*sn;
	err2*=err2;
	if(err2 == 0.0) err2=1.0;
		  Nprof++;
	X2_phase +=(line_mod[1][j]-flux+_par_[2][6])*(line_mod[1][j]-flux+_par_[2][6])/err2;
	//X2+= (line_mod[1][j]-flux+_par_[2][6])*(line_mod[1][j]-flux+_par_[2][6])/err2;
#ifdef VERBOSE_SLOW_ACTION
        cout << line_mod[0][j] << " " << line_mod[1][j] << " " 
             << _input_.dados[i].le1(j) << " " << flux << " " << sqrt(err2) <<endl;
#endif
      }
      else if(MLProf.le2(j) > 0) aborta("Problema no mapemento de pixels");
      
    }

#ifdef VERBOSE_SLOW
    cout<<endl;
#endif
	  X2 += X2_phase/Nprof;
	  NP+=Nprof; //double(_input_.dados[i].size());
	  //cout << X2_phase/Nprof << endl;
    //exit(0);
  }
  X2/=_input_.fase.size();
  //exit(0);
  if((_X2_ < 0) ||(_X2_ > X2)) _X2_=X2; //guarga o menor valor de X2 em _X2_
  _X2c_=X2;
#ifdef VERBOSE_SLOW
  cout << "# - X2_slow ended...\n";
#endif

  return X2;
}

//Slow functions apply instrumental blurring! Use it if instrumental bluring is considerable!
double DataKeeper::X2_slow_2(const gsl_vector &g)
{
	//Copia nova distribuicao de brilho para estrela...
	_estrela_.NewI(g);
#ifdef VERBOSE_SLOW
	cout << "# - X2_slow started...\n";
#endif
	double X2 = 0;
	double NP = 0;

	double X2_phase = 0.0;
	
	vector<double>& vfases = _input_.fase;
	
	double &sn=_par_[0][3];//Sinal-Ruido
	
	/* Mudo funcao que cuida de erro para nao aplicar blurringg em caso de problemas*/
	gsl_error_handler_t * old_handler = gsl_set_error_handler (&my_handler);


	for(int i = 0; i < vfases.size(); i++){
		// _input_.dados.col1() -> devolve referencia para vetor da coluna 1
		vector<double> lin_form(3);
		lin_form[0]= -_par_[1][4];
		lin_form[1]=  _par_[1][4];
		lin_form[2]=  _par_[1][5];
		Histograma MLProf = _estrela_.rad_prof(vfases[i], _input_.dados[i].col1(), _pixmap_[i]);
		
		// aplicando blurring instrumental
		vector<vector<double> > line_mod(2); //,MLProf.size());
		line_mod[0].resize(MLProf.size());
		line_mod[1].resize(MLProf.size());
		//cout << "#- Aqui 0\n";
		for(int ii = 0; ii < MLProf.size(); ii++){
			line_mod[0][ii] = MLProf.le1(ii);
			line_mod[1][ii] = MLProf.le2(ii);
		}
		double norm_fac = 0.0; //Fator de normalizacao de fluxo (Fdata/Fmod)
		//se tiver blurring aplica...
		if(_par_[1][6] > _par_[1][5]){
			Convolve ff(line_mod,_f2,line_mod[0].size(),_f2[0].size());
			ff.xmin(line_mod[0][0]);
			ff.xmax(line_mod[0][line_mod[0].size()-1]);
			norm_fac = _input_.fluxo[i]/_mlflux[i]; //ff.integrate1(); //_mlflux[i];
			
			//cout << i << " - aqui!\n";
			gsl_function F;
			F.function = &cnv_fnc;
			F.params = &ff;
			GSL_Integration itr_conv(F);
			itr_conv.xmin(ff.xmin());
			itr_conv.xmax(ff.xmax());
			//itr_conv.xmin(line_mod[0][4]); //(ff.xmin());
			//itr_conv.xmax(line_mod[0][line_mod[0].size()-4]); //(ff.xmax());
			
			for (int ii = 0; ii < line_mod[0].size(); ii++){
				ff.px(line_mod[0][ii]);
				line_mod[1][ii] = itr_conv.integrate_qag()/_par_[1][6]/sqrt(M_PI)*norm_fac;
				itr_conv.clear();
				if (global_my_gsl_err){
					cerr << "# - ERROR CALLED\n";					
//					string err = "lixo" ; 
//					write_res(err);
//					exit(1);
					
					break;
				}
				//cout << line_mod[0][ii]
				//<< " " << line_mod[1][ii] << endl;
			}
			ff.clean();
		}
		else{ //(global_my_gsl_err){
			//cerr << "# - SKIPING CONVOLUTION\n";
			global_my_gsl_err = false;
			Convolve ff(line_mod,_f2,line_mod[0].size(),_f2[0].size());
			ff.xmin(line_mod[0][0]);
			ff.xmax(line_mod[0][line_mod[0].size()-1]);
			norm_fac = _input_.fluxo[i]/ff.integrate1(); //_mlflux[i];
			
			for (int ii = 0; ii < line_mod[0].size(); ii++){
				line_mod[1][ii] *= _input_.fluxo[i]/ff.integrate1();
			}
			
		}
		
	

#ifdef VERBOSE_SLOW
			cout << "# - X2_slow - Convolution ended...\n";
#endif			
		
		//excluo os 2 primeiros e os 2 ultimos pontos do espectro devido a aplicacao do blurring!
#ifdef DEBUG_X2
		ostringstream fname;
		fname << "rprof_" << i << ".dat";
		ofstream filecprof(fname.str().c_str(),ios_base::app);
#endif
		for(int j = _input_.min_x[i]; j <= _input_.max_x[i]; j++){
			
			double flux = _input_.dados[i].le2(j);
			double err2 = _input_.err[i].le2(j)*sn;
#ifdef DEBUG_X2			
			filecprof << line_mod[0][j] << " "
			<< line_mod[1][j] << " " << flux+_par_[2][6] << endl;
#endif						
			X2_phase +=(line_mod[1][j]-flux+_par_[2][6])*(line_mod[1][j]-flux+_par_[2][6])/err2/err2;
			
		}
#ifdef DEBUG_X2
		filecprof.close();
#endif
		X2 += (X2_phase/(double)(_input_.max_x[i]-_input_.min_x[i])); 
		X2_phase = 0;
	}

	/* restore original handler */
	gsl_set_error_handler (old_handler);


	X2/=_input_.fase.size();
	//exit(0);
	if((_X2_ < 0) ||(_X2_ > X2)) _X2_=X2; //guarga o menor valor de X2 em _X2_
	_X2c_=X2;
#ifdef VERBOSE_SLOW
	cout << "# - X2_slow ended...\n";
#endif
	
	return X2;
}

/*************
X2_slow_3:
 Slow functions apply instrumental blurring! 
 Use it if instrumental bluring is considerable!
 This function combine profiles spaning the exposure time to account for orbital smearing.
*************/

double DataKeeper::X2_slow_3(const gsl_vector &g)
{
	//Copia nova distribuicao de brilho para estrela...
	_estrela_.NewI(g);
#ifdef VERBOSE_SLOW
	cout << "# - X2_slow started...\n";
#endif
	double X2 = 0;
	//double NP = 0;
	
	vector<double>& vfases = _input_.fase;

	double* X2_phase = new double [vfases.size()];

	for(int i = 0; i < vfases.size(); i++)
		X2_phase[i]=0.;
	
	double &sn=_par_[0][3];//Sinal-Ruido
	
	/* Mudo funcao que cuida de erro para nao aplicar blurringg em caso de problemas*/
	//gsl_error_handler_t * old_handler = gsl_set_error_handler (&my_handler);
	int BadProof = 0;// Flag to count bad profiles in case of intrumental blurr
	
	try {
		
//#pragma omp parallel for	
	for(int i = 0; i < vfases.size(); i++){
		// _input_.dados.col1() -> devolve referencia para vetor da coluna 1
		vector<double> lin_form(3);
		lin_form[0]= -_par_[1][4];
		lin_form[1]=  _par_[1][4];
		lin_form[2]=  _par_[1][5];
		Histograma* Vec_MLProf = new Histograma [_Np_];

		for (int nhist = 0; nhist < _Np_; ++nhist) {
			double dphi = _Np_ > 1 ? _par_[0][5]/(_Np_-1) * nhist - _par_[0][5]/2.0 : 0.0;

			Vec_MLProf[nhist] = _estrela_.rad_prof(vfases[i] + dphi , _input_.dados[i].col1(), _pixmap_[i+nhist*vfases.size()]);

		}

		vector<vector<double> > line_mod(2); //,MLProf.size());
		line_mod[0].resize(Vec_MLProf[0].size());
		line_mod[1].resize(Vec_MLProf[0].size());
		
		double norm_fac = 0.0; //_input_.fluxo[i]/_mlflux[i]; //Fator de normalizacao de fluxo (Fdata/Fmod)

		
		for(int ii = 0; ii < Vec_MLProf[0].size(); ii++){
			line_mod[0][ii] = Vec_MLProf[0].le1(ii);
			line_mod[1][ii] = 0.0; //MLProf.le2(ii);

			for (int kk = 0; kk < _Np_; ++kk) {

				line_mod[1][ii] += Vec_MLProf[kk].le2(ii);
			}
			line_mod[1][ii] /= double(_Np_);

			norm_fac+=line_mod[1][ii];

		}
/*********************************************************************************************************		  
		 for (int ii = 0; ii < line_mod[0].size(); ii++){
		 line_mod[1][ii] *= _input_.fluxo[i]/norm_fac;
		 
		 }
 *********************************************************************************************************/
		bool GoodProof = true;// Flag to avoid bad profiles in case of intrumental blurr
		
		if(_par_[1][6] > _par_[1][5]){
	
			Convolve ff(line_mod,_f2,line_mod[0].size(),_f2[0].size());
			ff.xmin(line_mod[0][0]);
			ff.xmax(line_mod[0][line_mod[0].size()-1]);
			norm_fac = _input_.fluxo[i]/_mlflux[i]; //ff.integrate1(); //_mlflux[i];
													//cout << i << " - aqui!\n";
			gsl_function F;
			F.function = &cnv_fnc;
			F.params = &ff;
			GSL_Integration itr_conv(F);
				//itr_conv.xmin(ff.xmin());
				//itr_conv.xmax(ff.xmax());
            itr_conv.xmin(line_mod[0][4]); //(ff.xmin());
            itr_conv.xmax(line_mod[0][line_mod[0].size()-4]); //(ff.xmax());
			norm_fac = 0;
			double avg_itrerr = 0.;
			for (int ii = 4; ii < line_mod[0].size()-4; ii++){
				ff.px(line_mod[0][ii]);
				line_mod[1][ii] = itr_conv.integrate_cquad()/sqrt(M_PI)/_par_[1][6];
				norm_fac += line_mod[1][ii];
				avg_itrerr += itr_conv.integrationErr(); // *line_mod[1][ii]; // *norm_fac;
				itr_conv.clear();
					//cout << line_mod[0][ii]
					//<< " " << line_mod[1][ii] << endl;
			}
			if (avg_itrerr/line_mod[0].size() < 1.0){
				++BadProof;
				GoodProof = false;
			}
			ff.clean();
		}
		
		if (BadProof > vfases.size()/10){
			cout << "\nToo Many bad profile (" << BadProof << "/" << vfases.size() << ")" <<endl;
			string root = "/tmp/badprof.dat";
			write_res(root);
			throw "Bad Blurr";
		}
			//norm_fac /= double(Vec_MLProf[0].size());
		
		delete [] Vec_MLProf;
		
		if (GoodProof){

			norm_fac /= double(line_mod[0].size());
		
			for (int ii = 0; ii < line_mod[0].size(); ii++)
				line_mod[1][ii] *= _input_.fluxo[i]/_mlflux[i]; // norm_fac;
		
#ifdef VERBOSE_SLOW
		cout << "# - X2_slow - Convolution ended...\n";
#endif			
		
		//excluo os 2 primeiros e os 2 ultimos pontos do espectro devido a aplicacao do blurring!
#ifdef DEBUG_X2
		ostringstream fname;
		fname << "rprof_" << i << ".dat";
		ofstream filecprof(fname.str().c_str(),ios_base::app);
#endif
			for(int j = _input_.min_x[i]; j <= _input_.max_x[i]; j++){
			
				double flux = _input_.dados[i].le2(j);
				double err2 = _input_.err[i].le2(j)*sn;
#ifdef DEBUG_X2				
			filecprof << line_mod[0][j] << " "
			<< line_mod[1][j] << " " << flux+_par_[2][6] << endl;
#endif						
				if (line_mod[1][j] < __ZERO_MOD_FLUX__)
					X2_phase[i] +=(line_mod[1][j]-flux+_par_[2][6])*(line_mod[1][j]-flux+_par_[2][6])/err2/err2;

			}
#ifdef DEBUG_X2
		filecprof.close();
#endif
			X2_phase[i] /= (double)(_input_.max_x[i]-_input_.min_x[i]);
		//X2 += (X2_phase/(double)(_input_.max_x[i]-_input_.min_x[i])); 
		//X2_phase = 0;
		}

	}
		
	} catch (char* code) {
		cerr << "Too much bad profiles produced due to blurring! You should try turning blurr off.\n";
		exit(-1);
	}

	/* restore original handler */
	//gsl_set_error_handler (old_handler);

	X2 = 0.;
	for(int i = 0; i < vfases.size(); i++)
		X2+=X2_phase[i];

	delete X2_phase;

	X2/=(double)_input_.fase.size();
	//exit(0);
	if((_X2_ < 0) ||(_X2_ > X2) && !isnan(X2)) _X2_=X2; //guarga o menor valor de X2 em _X2_
	_X2c_=X2;
#ifdef VERBOSE_SLOW
	cout << "# - X2_slow ended...\n";
#endif
	
	return X2;
}


void DataKeeper::dX2_slow(const gsl_vector &x, gsl_vector &rdX2)
{
	/* Mudo funcao que cuida de erro para nao aplicar blurringg em caso de problemas*/
	//gsl_error_handler_t * old_handler = gsl_set_error_handler (&my_handler);

  _estrela_.NewI(x);
#ifdef VERBOSE_SLOW
  cout << "# - dX2_slow started...\n";
#endif

  //gsl_vector* rdX2 = gsl_vector_alloc(x.size); //vetor de retorno
  
  double &sn=_par_[0][3];//Sinal-Ruido
  double ABS_FLAG = _par_[0][4]; //Flag de absorcao
  //float ABS_FLAG = 1.0; // -1 p/ emissao e 1 para absorcao
	if (abs(ABS_FLAG) != 1.0) ABS_FLAG = ABS_FLAG/abs(ABS_FLAG);

  //Modeled Line Profile:
  vector<vector<double> > line_mod(_input_.fase.size());

  for(int k = 0; k < _input_.fase.size(); k++){
    Histograma MLProf = _estrela_.rad_prof(_input_.fase[k], _input_.dados[k].col1(), _pixmap_[k]);
    // aplicando blurring instrumental
    vector<double> velvec;
    velvec.resize(MLProf.size());
    line_mod[k].resize(MLProf.size());
    //cout << "#- Aqui 0\n";
    vector<vector<double> > line_mod2(2);
    line_mod2[0].resize(MLProf.size());
    line_mod2[1].resize(MLProf.size());

    for(int ii = 0; ii < MLProf.size(); ii++){
      velvec[ii] = MLProf.le1(ii);
      line_mod[k][ii] = MLProf.le2(ii);
      line_mod2[0][ii]=velvec[ii];
      line_mod2[1][ii]=line_mod[k][ii];
    }
    //cout << "#- Aqui 1\n";
    //tmp_vec=vector<double> (lp.col2());
    //line_mod[1].push_back(tmp_vec);
    //se tiver blurring aplica...
    double norm_fac;
	  
	  	  
    if(_par_[1][6] > _par_[1][5]){
      Convolve ff(line_mod2,_f2,line_mod[k].size(),_f2[0].size());
      ff.xmin(velvec[0]);
      ff.xmax(velvec[velvec.size()-1]);
      norm_fac= _input_.fluxo[k]/ff.integrate1(); //_mlflux[k]; /
      gsl_function F;
      F.function = &cnv_fnc;
	  F.params = &ff;
      GSL_Integration itr_conv(F);
	  // itr_conv.xmin(_input_.dados[k].le1(0)); //(ff.xmin());
      // itr_conv.xmax(_input_.dados[k].le1(_input_.dados[k].size()-1)); //(ff.xmax());
      itr_conv.xmin(_input_.dados[k].le1(_input_.min_x[k])); //(line_mod[0][4]); //(ff.xmin());
      itr_conv.xmax(_input_.dados[k].le1(_input_.max_x[k])); //(line_mod[0][line_mod[0].size()-4]); //(ff.xmax());
		
		
	  for (int ii = 0; ii < velvec.size(); ii++){
        ff.px(velvec[ii]);
        line_mod[k][ii] = itr_conv.integrate_qag()/_par_[1][6]/sqrt(M_PI)*norm_fac;
        itr_conv.clear();
		  if (global_my_gsl_err){
			  cerr << "# - ERROR CALLED\n";
			  break;
		  }
		  
        //if(ii == velvec.size()-1)itr_conv.clean();
        //cout << line_mod[0][i]
        //<< " " << line_mod[1][i] << " " 
        //<< itr_conv.integrate_qag()/2.0 << endl;
      }
      ff.clean();
    }
		else { //if (_par_[1][5] > 0 || global_my_gsl_err){
		  //cerr << "# - dX2 SKIPING CONVOLUTION\n";
		  global_my_gsl_err = false;
		 
		  
		Convolve ff(line_mod2,_f2,line_mod[k].size(),_f2[0].size());
		ff.xmin(velvec[0]);
		ff.xmax(velvec[line_mod[0].size()-1]);
		norm_fac = _input_.fluxo[k]/ff.integrate1(); //_mlflux[i];
		
		for (int ii = 0; ii < velvec.size(); ii++){
			line_mod[k][ii] *= _input_.fluxo[k]/ff.integrate1();
		}
		
	}
	  
	  
  }

  for(int i = 0; i < x.size; i++){

    double X2 = 0;
    for(int k = 0; k < _input_.fase.size(); k++){
      double vis = _estrela_.VFunc(i, _input_.fase[k]);
      if(vis > 0){ 
	// Posicao do pixel no histograma de velocidade 
	int HistPixPos = _pixmap_[k][i];
//excluo os 2 primeiros e os 2 ultimos pontos do espectro devido a aplicacao do blurring!
	//if((2 <= HistPixPos) && (HistPixPos < _input_.dados.size()-2)){
	  double flux = _input_.dados[k].le2(HistPixPos);
	  double err2 = _input_.err[k].le2(HistPixPos)*sn; //flux/sn;
	  //double err2 = _input_.dados[k].le2(HistPixPos)*sn; //flux/sn;
	  //err2*=err2;
	  if(err2 == 0.0) err2=1.0;
	  X2+= 2.0*vis*( line_mod[k][HistPixPos] - flux - _par_[2][6] ) / err2;
	 
	//}	
      }
    }
	  //cout << X2 << endl;
	  //Para absorcao tem que t

	  gsl_vector_set(&rdX2, i, (ABS_FLAG*X2/double(_input_.fase.size())));

  }
	
#ifdef VERBOSE_SLOW
  cout << "# - dX2_slow ended...\n";
#endif

  //return rdX2;
}



void DataKeeper::dX2_slow_3(const gsl_vector &x, gsl_vector &rdX2)
{
	
	_estrela_.NewI(x);
#ifdef VERBOSE_SLOW
	cout << "# - dX2_slow started...\n";
#endif
	
	
	double &sn=_par_[0][3];//Sinal-Ruido
	double ABS_FLAG = _par_[0][4]; //Flag de absorcao

	if (abs(ABS_FLAG) != 1.0) ABS_FLAG = ABS_FLAG/abs(ABS_FLAG);
	
	//Modeled Line Profile:
	vector<vector<double> > line_mod(_input_.fase.size());
	int BadProof = 0;// Flag to count bad profiles in case of intrumental blurr	
	bool* GoodProof = new bool [_input_.fase.size()];// Flag to avoid bad profiles in case of intrumental blurr	
#pragma omp parallel for
	for(int k = 0; k < _input_.fase.size(); k++){
		Histograma* MLProf = new Histograma [_Np_];
		
		for (int nhist = 0; nhist < _Np_; ++nhist) {
			double dphi = _Np_ > 1 ? _par_[0][5]/(_Np_-1) * nhist - _par_[0][5]/2.0 : 0;
			
			MLProf[nhist] = _estrela_.rad_prof(_input_.fase[k] + dphi , _input_.dados[k].col1(), _pixmap_[k+nhist*_input_.fase.size()]);
		}

		
		// aplicando blurring instrumental
		vector<double> velvec;
		
		velvec.resize(MLProf[0].size());
		
		line_mod[k].resize(MLProf[0].size());
		

		vector<vector<double> > line_mod2(2);
		line_mod2[0].resize(MLProf[0].size());
		line_mod2[1].resize(MLProf[0].size());
		
		double norm_fac = 0.0;
		

		for(int ii = 0; ii < MLProf[0].size(); ii++){
			velvec[ii] = MLProf[0].le1(ii);
			line_mod[k][ii] = 0.0;
			for (int nprof = 0; nprof < _Np_;++nprof)
				line_mod[k][ii] += MLProf[nprof].le2(ii);
			
			line_mod[k][ii] /= double(_Np_);
			
			norm_fac+=line_mod[k][ii];
			
			line_mod2[0][ii]=velvec[ii];
			line_mod2[1][ii]=line_mod[k][ii];
		}
		
		delete [] MLProf;
/**********************************************************************************************************/
		
		GoodProof[k] = true;
		
		if(_par_[1][6] > _par_[1][5]){

			Convolve ff(line_mod2,_f2,line_mod2[0].size(),_f2[0].size());
			ff.xmin(line_mod2[0][0]);
			ff.xmax(line_mod2[0][line_mod[0].size()-1]);
			//norm_fac = _input_.fluxo[k]/_mlflux[k]; //ff.integrate1(); //_mlflux[i];
													//cout << i << " - aqui!\n";
			gsl_function F;
			F.function = &cnv_fnc;
			F.params = &ff;
			GSL_Integration itr_conv(F);
				//itr_conv.xmin(ff.xmin());
				//itr_conv.xmax(ff.xmax());
            itr_conv.xmin(line_mod2[0][4]); //(ff.xmin());
            itr_conv.xmax(line_mod2[0][line_mod2[0].size()-4]); //(ff.xmax());

			double avg_itrerr = 0.;
			for (int ii = 4; ii < line_mod[0].size()-4; ii++){
				ff.px(line_mod2[0][ii]);
				line_mod2[1][ii] = itr_conv.integrate_cquad()/sqrt(M_PI)/_par_[1][6];
				//norm_fac += line_mod2[1][ii];
				avg_itrerr += itr_conv.integrationErr(); // *line_mod[1][ii]; // *norm_fac;
				itr_conv.clear();
				
					//cout << line_mod[0][ii]
					//<< " " << line_mod[1][ii] << endl;
			}

			if (avg_itrerr/line_mod2[0].size() < 1.0){
				++BadProof;
				GoodProof[k] = false;
			}
			ff.clean();
		

        }

/**********************************************************************************************************/
		if (GoodProof[k]){

			norm_fac /= double(line_mod2[0].size());

			for(int ii = 0; ii < line_mod2[0].size(); ii++)
				line_mod[k][ii] = line_mod2[1][ii]/_input_.fluxo[k]/_mlflux[k];		


		}
	}
	
	for(int npix = 0; npix < x.size; ++npix){ //iterando em pixel
		
		double dX2 = 0.0; // gradiente do pixel
		
		for(int npha = 0; npha < _input_.fase.size(); ++npha){ //iterando em fase
			
			if (_pixmap_[npha+int(_Np_/2)*_input_.fase.size()][npix] > 0 && GoodProof[npha]) { 
				// se pixel for visivel no perfil do meio eu considero, senao pulo!
				
				for (int nhist = 0; nhist < _Np_; ++nhist) {
					double dphi = _Np_ > 1 ? _par_[0][5]/(_Np_-1) * nhist - _par_[0][5]/2.0 : 0;
					double vis = _estrela_.VFunc(npix, _input_.fase[npha]+dphi);
					double dflux = _input_.dados[npha].le2(_pixmap_[npha+nhist*_input_.fase.size()][npix]);
					double derro = _input_.err[npha].le2(_pixmap_[npha+nhist*_input_.fase.size()][npix]);
					double& mflux = line_mod[npha][_pixmap_[npha+nhist*_input_.fase.size()][npix]];
					dX2 += vis > 0.0 ? 2.0 * ( mflux - dflux ) / derro / sn : 0.0;
				}
			}
			
/************************************************************************
			double vis = _estrela_.VFunc(i, _input_.fase[k]);
			
			if(vis > 0){ 
				// Posicao do pixel no histograma de velocidade 
				int HistPixPos = _pixmap_[k][i];
				//excluo os 2 primeiros e os 2 ultimos pontos do espectro devido a aplicacao do blurring!
				//if((2 <= HistPixPos) && (HistPixPos < _input_.dados.size()-2)){
				double flux = _input_.dados[k].le2(HistPixPos);
				double err2 = _input_.err[k].le2(HistPixPos)*sn; //flux/sn;
				//double err2 = _input_.dados[k].le2(HistPixPos)*sn; //flux/sn;
				//err2*=err2;
				if(err2 == 0.0) err2=1.0;
				X2+= 2.0*vis*( line_mod[k][HistPixPos] - flux - _par_[2][6] ) / err2;
				
				//}	
			}
 ************************************************************************/
		}
		//cout << X2 << endl;
		//Para absorcao tem que t
		
		gsl_vector_set(&rdX2, npix, (ABS_FLAG*dX2/double(_input_.fase.size())));
		
	}
	
	
#ifdef VERBOSE_SLOW
	cout << "# - dX2_slow ended...\n";
#endif
	
	//return rdX2;
}


/*****************
 
 END OF CHI-SQUARE DEFINITIONS
 
 *****************/

void DataKeeper::dS(const gsl_vector &x,gsl_vector &rdS)
{
	double s_ = 0.0;
	double soma = 0.0;
	
	for (int k = 0; k < x.size; k++) {
		double gx = gsl_vector_get(&x, k);
		s_ -= gx * gx * log(gx * gx);
		soma += ( gx * gx );
	}
	s_ /= soma;
	
	for (int k = 0; k < x.size; k++) {
		gsl_vector_set( &rdS, k, log( gsl_vector_get(&x, k) * gsl_vector_get(&x, k) ) + s_ / -soma);
	}
}
/*
{
	//double DX = 5e-8
  //_estrela_.NewI(x);
  mval(x);
  //gsl_vector* rdS = gsl_vector_alloc(x.size);
	const double min = 0.;
	const double max = 2.;
	
	//cout << "# - Aqui S(0)\n";
	vector<double> img (x.size,0.);
	//vector<double> dimg (x.size,0.);
	
	for (int i = 0; i < x.size; i++) {
		img[i] = gsl_vector_get(&x, i);
		//img[i] = gsl_vector_get(&x, i);
	}
	
	vector<double> fmt(3);
	fmt[0] = min; // gsl_vector_min(&g);
	fmt[1] = ( gsl_vector_max(&x) < max ) ? max : gsl_vector_max(&x);
	fmt[2] =  (fmt[1] - fmt[0]) / sqrt(double(x.size));
	
	Histograma p_hist(img,fmt);
	
	for(int i = 0; i < x.size; i++){
		int iindex = -1;
		//cout << gsl_vector_get(&x,i) <<' ' << mval() << ' '<<-log(gsl_vector_get(&x,i)/mval())<<endl;
		//gsl_vector_set(&rdS,i,log(gsl_vector_get(&x,i)/mval()));
		
		for (int k = 0; k < p_hist.size(); k++) {
			if ( ( p_hist.le1(k)-fmt[2]/2. <= gsl_vector_get(&x,i) ) && ( gsl_vector_get(&x,i) < p_hist.le1(k)+fmt[2]/2. ) ) {
				iindex = k;
				break;
			}
		}
		if ( ( iindex > 0 ) && ( p_hist.le2(iindex) > 0) ) {
			//gsl_vector_set(&rdS,i,gsl_vector_get(&x,i)*log(1.0/p_hist.le2(iindex)));
			//double dS_ = 1./(gsl_vector_get(&x,i)*p_hist.le2(iindex))*
			//(1.+log10(1./p_hist.le2(iindex))-1./p_hist.le2(iindex)-1./p_hist.le2(iindex)*log10(1./p_hist.le2(iindex)));
			gsl_vector_set(&rdS,i,log(1.0/p_hist.le2(iindex)));
			//cout << iindex  << " " << p_hist.le1(iindex) << " " << p_hist.le2(iindex) << " " << log10(1.0/p_hist.le2(iindex)) << " " << endl ;
		}
		else {
			cout << "# - WARNING: Index out of range ...\n";
			cout << "# - WARNING: I = " << gsl_vector_get(&x, i) << endl;
			cout << "# - WARNING: k = " << iindex << "[" << p_hist.size() << "] -> " << p_hist.le1(iindex)-fmt[2]/2. << " <= I < " << p_hist.le1(iindex)+fmt[2]/2. << endl<<endl;
		}


		//gsl_vector_set(&rdS,i, 2.0*( gsl_vector_get(&x,i)-mval() ) );
		// o sinal é negativo para que, qdo o pixel tiver mais riblho que o valor medio ele tenha que
		// ser diminuido e se ele for menor ele tem que crescer...
	}
	//exit(0);
  //return rdS;
}
*/

void DataKeeper::dS_gauss_fast(const gsl_vector &x,gsl_vector &rdS)
{
	//_estrela_.NewI(g);
	//double S = 0.0;
	//mval(x);
	/*
	for(int i = 0; i < x.size; i++){
		double Mj = gsl_vector_get(_Default_Map,i);
		double Dj = gsl_vector_get(&x,i); //Dj = mval();
		gsl_vector_set(&rdS,i, 2.0*( Dj-Mj ) );
	}
	*/
	
	double s_ = 0.0;
	double soma = 0.0;
	double somadef = 0.0;
	for (int k = 0; k < x.size; k++) {
		double gx = gsl_vector_get(&x, k);
		double gdef = gsl_vector_get(_Default_Map,k);
		s_ -= gx * gx * log( (gx * gx) / (gdef * gdef) );
		soma += ( gx * gx );
		somadef += ( gdef * gdef );
	}
	s_ /= soma;
	s_ += log( soma / somadef );
	
	for (int k = 0; k < x.size; k++) {
		gsl_vector_set( &rdS, k, log( (gsl_vector_get(&x, k) * gsl_vector_get(&x, k)) / ( gsl_vector_get(_Default_Map,k) * gsl_vector_get(_Default_Map,k) ) ) + s_ / -soma);
	}
	
}


void DataKeeper::show_par(ostream& out)
{
  vector<vector<string> > par_def (4);

  const int p = 0;
  const int p0 = 1;
  const int p1 = 2;
  const int p2 = 3;

  par_def[p].push_back("SIZE = ");
  par_def[p].push_back("Xaim = ");
  par_def[p].push_back("FWHM = ");
  par_def[p].push_back("S/N = ");
  par_def[p].push_back("ABS_FLAG = ");
  par_def[p].push_back("texp/Porb = ");

  par_def[p0].push_back("i = ");
  par_def[p0].push_back("M1 = ");
  par_def[p0].push_back("M2 = ");
  par_def[p0].push_back("P = ");
  //par_def[p0].push_back("S/N = ");
  par_def[p0].push_back("Vmax = ");
  par_def[p0].push_back("lres = ");
  par_def[p0].push_back("lres2 = ");

  par_def[p1].push_back("f2 = ");
  par_def[p1].push_back("grav_coef = ");
  par_def[p1].push_back("limb_law = ");
  par_def[p1].push_back("nspot = ");
  par_def[p1].push_back("T_pole = ");
  par_def[p1].push_back("ffac = ");
  par_def[p1].push_back("fd = ");

  for(int i = 0; i < _par_[p2].size(); i++){
    ostringstream str0;
    str0 << "x" << i << " = ";
    par_def[p2].push_back(str0.str());
  }

  out << "# Tiago Ribeiro - LinuxC++ - 2006\n#";

  for(int i = 0; i < _par_.size(); i++){
    int nc = 0;
    if(_par_[i].size() != par_def[i].size()){
      cout << "WARNING - Definicoes diferentes entre par e par_def.\n"
	  << "par = " << i << " "  << _par_[i].size() << endl
	  << "par_def = " << i << " "  << par_def[i].size() << endl;
      if(_par_[i].size() > par_def[i].size()) nc = par_def[i].size();
      else nc = _par_[i].size();
    }
    else nc = _par_[i].size();
    for(int j = 0; j < nc; j++){
      if( i < 4 ) out << par_def[i][j];
      else out << par_def[3][j];
      
      if( i == 1 && j == 0) out << _par_[i][j]*180/M_PI << ' ';
      else out << _par_[i][j] << ' ';
    }
    out << "\n#";
  }
//  out << "V(Km/s)  Flux\n";

}

double DataKeeper::S(const gsl_vector &g)
{
	double s_ = 0.0;
	double soma = 0.0;
	double somadef = 0.0;
	mval(g);
	double gdef = mval();	
	for (int k = 0; k < g.size; k++) {
		double gx = gsl_vector_get(&g, k);

		
		//cout << gx * gx * log( (gx * gx) / (gdef * gdef) )  << " " << gx << " " << gdef << endl;
		
		s_ -= gx * gx * log( (gx * gx) / (gdef * gdef) );
		soma += ( gx * gx );
		somadef += ( gdef * gdef );
	}
	s_ /= soma;
	s_ += log( soma / somadef );
	
	if((_S_ <= 0) || (_S_ > s_) && !isnan(s_))_S_ = s_;
	
	_Sc_=s_;
	if(isnan(s_)) s_ = 0;
	//cout << "\n## S[INTERNAL] = " << s_ << endl;;
	return -sqrt(s_*s_);
	
}
/*
{
	const double min = 0.;
	const double max = 2.;
	
	//cout << "# - Aqui S(0)\n";
	vector<double> img (g.size,0.);
	
	for (int i = 0; i < g.size; i++) {
		img[i] = gsl_vector_get(&g, i);
	}
	
	vector<double> fmt(3);
	fmt[0] = min; // gsl_vector_min(&g);
	fmt[1] = ( gsl_vector_max(&g) < max ) ? max : gsl_vector_max(&g);
	fmt[2] =  (fmt[1] - fmt[0]) / sqrt(double(g.size));
	
	Histograma p_hist(img,fmt);

	double s_ = 0.; // stores entropy
	double ns_ = 0.;
	
	for (int i = 0; i < p_hist.size(); i++) {
		if (p_hist.le2(i) > 0.0) {
			s_ += 1./p_hist.le2(i) * log10(1./p_hist.le2(i));
			ns_++;
		}
		//cout <<	p_hist.le1(i) << " " << p_hist.le2(i) << endl;
	}
	//exit(0);
	(ns_ > 0.) ? s_ /= ns_ : s_ = s_ ;
	
	if((_S_ <= 0) || (_S_ > s_))_S_ = s_;

	_Sc_=s_;

	return s_;
}

/*
 double DataKeeper::S_OLD(const gsl_vector &g)
{
  //_estrela_.NewI(g);
  double S = 0;
  mval(g);
  
  for(int i = 0; i < g.size; i++){
    double Mj = gsl_vector_get(&g,i);
    double Dj = mval();
    double logMjDj = log10(Mj/Dj); //(Mj > Dj) ? log(Mj/Dj) : 0;

    //S += Mj-Dj-Mj*logMjDj;
	  S -= (Mj-Dj)*(Mj-Dj); //-abs(Mj/Dj * logMjDj);
	  //if (Mj-Dj-Mj*logMjDj 
	  //S+=logMjDj;
    //cout << log(Mj/Dj) << " " << Mj << " " << Dj << endl;
  }
  //cout << " "  << S << endl;
  //exit(0);
  S/=double(g.size);
  if((_S_ <= 0) || (_S_ > S))_S_ = S;
  _Sc_=S;
  return S;
}
*/
 
double DataKeeper::S_gauss_slow(const gsl_vector &g,const double FWHM)
{
	//_estrela_.NewI(g);
	double S = 0.0;
	mval(g);
	surface Srf = image();
	for(int i = 0; i < g.size; i++){
		double* pixlist;
		//cout << "# aqui 1\n";
		int npix = _estrela_.close_pix(Srf.px[i],Srf.py[i],Srf.pz[i],FWHM, pixlist);
		//cout << "# aqui 2\n";
		double Mj = 0.0;
		for( int j = 0; j < npix; j++)
			Mj += gsl_vector_get(&g,pixlist[j]);
		delete pixlist;
		//cout << "# aqui 3 - " << i << " / " << g.size << " / " << npix << "\n";
		Mj /= double(npix);
		double Dj = gsl_vector_get(&g,i); //Dj = mval();
		//double logMjDj = log10(Mj/Dj); //(Mj > Dj) ? log(Mj/Dj) : 0;
		
		//S += Mj-Dj-Mj*logMjDj;
		S -= (Mj-Dj)*(Mj-Dj); //-abs(Mj/Dj * logMjDj);
		//if (Mj-Dj-Mj*logMjDj 
		//S+=logMjDj;
		//cout << log(Mj/Dj) << " " << Mj << " " << Dj << endl;
	}
	//cout << " "  << S << endl;
	//exit(0);
	S/=double(g.size);
	if((_S_ <= 0) || (_S_ > S))_S_ = S;
	_Sc_=S;
	return S;
}

double DataKeeper::S_gauss_fast(const gsl_vector &g)
{
	double s_ = 0.0;
	double soma = 0.0;
	double somadef = 0.0;
	for (int k = 0; k < g.size; k++) {
		double gx = gsl_vector_get(&g, k);
		double gdef = gsl_vector_get(_Default_Map,k);
		
		//cout << gx * gx * log( (gx * gx) / (gdef * gdef) )  << " " << gx << " " << gdef << endl;
		
		s_ -= gx * gx * log( (gx * gx) / (gdef * gdef) );
		soma += ( gx * gx );
		somadef += ( gdef * gdef );
	}
	s_ /= soma;
	s_ += log( soma / somadef );
	
	if((_S_ <= 0) || (_S_ > s_) && !isnan(s_))_S_ = s_;
	
	_Sc_=s_;
	if(isnan(s_)) s_ = 0;
	//cout << "\n## S[INTERNAL] = " << s_ << endl;;
	return -sqrt(s_*s_);
	
	/*
	 const double min = 0.;
	 const double max = 2.;
	 
	 //cout << "# - Aqui S(0)\n";
	 vector<double> img (g.size,0.);
	 vector<double> dimg (g.size,0.);
	 
	 for (int i = 0; i < g.size; i++) {
	 img[i] = gsl_vector_get(&g, i);
	 dimg[i] = gsl_vector_get(_Default_Map,i);
	 }
	 
	 vector<double> fmt(3);
	 fmt[0] = min; // gsl_vector_min(&g);
	 fmt[1] = ( gsl_vector_max(&g) < max ) ? max : gsl_vector_max(&g);
	 fmt[2] =  (fmt[1] - fmt[0]) / sqrt(double(g.size));
	 
	 Histograma p_hist(img,fmt);
	 Histograma dp_hist(dimg,fmt);
	 
	 double s_ = 0.; // stores entropy
	 
	 for (int i = 0; i < p_hist.size(); i++) {
	 if ( (p_hist.le2(i) > 0) && (dp_hist.le2(i) > 0)) {
	 s_ += 1./p_hist.le2(i) * log10(dp_hist.le2(i) / p_hist.le2(i));
	 }
	 //cout <<	p_hist.le1(i) << " " << p_hist.le2(i) << endl;
	 }
	 //exit(0);
	 return s_;
	 
	/*	
	 //_estrela_.NewI(g);
	 double S = 0.0;
	 mval(g);
	 for(int i = 0; i < g.size; i++){
	 double Mj = gsl_vector_get(_Default_Map,i);
	 double Dj = gsl_vector_get(&g,i); //Dj = mval();
	 S -= (Mj-Dj)*(Mj-Dj); //-abs(Mj/Dj * logMjDj);
	 }
	 S/=double(g.size);
	 if((_S_ <= 0) || (_S_ > S))_S_ = S;
	 _Sc_=S;
	 return S;
	 
*/
 }


vector<vector<double> > org_par(vector<vector<double> > &in_par)
{

  vector<vector<double> > saida;

  // Infelismente nao existe maneira automatica de fazer isso. :(

  //  cout << endl;

  // Neste caso so preciso "deslocar" o vetor in_par de um e ja estou com o 
  // formato desejado!

  saida.push_back(in_par[1]); 
  saida.push_back(in_par[2]); 
  saida.push_back(in_par[3]); 
  //cout << saida[0][0] << endl;
  return saida;
}


int definicoes_obscurecimento(const int law)
{
  vector<vector<double> > par;

  Light_Curve teste(par);
  
  return teste.npar_obs_borda(law);
}


